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Abstract. It is quite common that several different phases exist simuhaneously in 
a system of trapped quantum gases of uhra-cold atoms. One such example is the 
strongly-interacting Fermi gas with two imbalanced spin species, which has received a 
great amount of attention due to the possible presence of exotic superfluid phases. By 
employing novel numerical techniques and algorithms, we self-consistently solve the 
Bogoliubov de- Gennes equations, which describe Fermi superfluids in the mean- field 
framework. From this study, we investigate the novel phases of spin-imbalanced Fermi 
gases and examine the validity of the local density approximation (LDA), which is often 
invoked in the extraction of bulk properties from experimental measurements within 
trapped systems. We show how the validity of the LDA is affected by the trapping 
geometry, number of atoms and spin imbalance. 
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1. Introduction 

Interest in spin-imbalanced Fermi superfluids dates back over a half century to 
the appearance of the Bardeen-Cooper-Schrieffer (BCS) theory. Several pioneers 
theoretically considered the fate of a superconductor in the presence of magnetic 
impurities which could disrupt the spin balance [H [21 |3l H] prescribed by BCS. In 
the ensuing decades the BCS picture has emerged as an important paradigm in many 
branches of physics [5], such as nuclear physics, quantum chromo dynamics, and ultra- 
cold atoms. The ultra-cold atoms, in particular, provide a highly controllable clean 
system to investigate the effect of spin imbalance on the nature of fermionic pairing. 
Over the past few years, experiments on polarized Fermi gases at Rice University [SllTllH], 
MIT P [ini ttH 112] and ENS p] have stimulated a flurry of theoretical activity. These 
efforts represent an important area within the general goal of emulating complicated 
many-body systems using cold atoms. 

All cold atom experiments are necessarily performed in the presence of trapping 
potentials that hold the atoms together in an inhomogeneous environment. In order to 
extract the bulk properties of the system (e.g, the equations of state) from measurement 
on a trapped sample, a local density approximation (LDA) is employed to account for 
the effects of the trapping potential. The LDA states that the system can be treated 
locally as a part of an infinite system. In other words, on a length scale which is short 
in comparison to the spatial variation of the external potential, the external potential 
Kxt(r) essentially acts as an offset to the chemical potential. In practice, one defines a 
local chemical potential yu(r) = /i — Vcxt{i^), with fi being the global chemical potential. 
In the context of cold atoms, the LDA is often an accurate approximation when there 
are no phase boundaries present. However, because of the large change in the density 
which occurs from the center to the edge of the trap, there is often more than one phase 
present and care needs to be taken in application of the LDA. This is especially true at 
phase boundaries where the LDA usually fails. Even in the circumstances where there 
is only one phase present in the trap, corrections to the LDA become crucial when the 
number of particles in the sample is small and finite-size effects are significant [H] . 

The trapped polarized Fermi gas represents just such a system — in all experimental 
observations, phase separation between a region with vanishing spin population and one 
with finite spin population has been observed. The purpose of this paper is to examine 
the validity of the LDA as the particle number and aspect ratio of the trap is changed. To 
this end we employ the Bogoliubov-de Gennes (BdG) equations |l5], which is a powerful 
mean-field tool particularly suitable for inhomogeneous Fermi superfluids, and has been 
recently adopted by many to study trapped ultra-cold Fermi gases [HI [T71 [HI [191 120] . 

In a recent paper [21] we analyzed discrepancies observed in pioneering experiments 
on polarized fermionic superfluids [9l [TOl 111 [Z] which appear to center on physics 
arising when the containing trap is highly elongated. In our analysis we found that 
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as the trap becomes increasingly elongated, the solutions of the BdG show a tendency 
towards metastable behavior which could lead to the observation of states which are 
not necessarily the ground state. However we also observed that one class of solutions 
with structure similar to the LDA solution was consistently the lowest in energy within 
our analysis. Our conclusions have since been confirmed by similar calculations [22] 
using a density functional formulation |23] which is more sophisticated than the BdG 
and accounts for quantum fluctuations and interactions within the normal fluid. In this 
paper we focus on this class of solutions and examine how well such solutions match 
those obtained from the LDA approximation as a function of particle number and trap 
aspect ratio. 

The content of the paper is organized as follows. In Sec. [2l we present the BdG 
formulation and describe the numerical techniques used to solve the BdG equations. In 
Sec. [3] we discuss our results for a relatively small number of particles and for different 
trapping geometry. In Sec. IH we focus on a very elongated cigar-like trap but vary the 
number of atoms. Finally, a concluding remark is presented in Sec. El 



2. The Bogoliubov De-Gennes treatment 

2.1. Formulation 

We consider a gas of spin-polarized fermionic atoms interacting through a contact 
potential {g^i^j S^i'f^i — ^j)) and confined to a harmonic trap defined in cylindrical 
coordinates {p,(f),z) by Vcxt{p,z) = y(a;^p^ + u'^z'^) with axial and radial frequencies 
denoted by {uz,uj±). We work at unitarity (a^ — )■ oo) and within a cigar-shaped 
trap with aspect ratio a = u±/ujz. This system of N = + atoms interacting 
through contact interaction is described by a Hamiltonian H = J dR {Hq + Hj) with 
non-interacting Hq and interacting Hj portions given by: 

Ho{R) = J2 + (P, z) - ^i^)^P^ , 

Hj{R) = -g^\{R)^l{R)^^{R)^^{R) , (1) 

where i}j<t{R) represents the fermionic field operators, m the mass and /lo- the chemical 
potential of atomic species with spin a. The coupling constant is defined as g = 
Henceforth, we work in trap units for which: m = Uz = h = ks = 1- This implies that 
energies will be measured in units of hiOz, lengths in units of Iq = \J^^ and temperature 
(T) in units of hUz/kB ■ 

The Hamiltonian ([T]) will be treated within the mean-field BdG approximation for 
which there are many excellent references [151 1211 [25]. Here we simply state the BdG 
equations for the pair wave functions Uj{R) and Vj{R) which decouple H: 





Uj 




Uj 











(2) 
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In the above coupled set of equations, Uj{R) and Vj{R) are two components of the quasi- 
particle wavefunction associated with energy Ej. The single particle Hamiltonian is 
defined in our trap units by: 

K(^) = -VV2 + Kxt + (?Pa-/i<x (3) 

and includes the trapping potential, the chemical potential fi^^ and the Hartree mean- 
field potential is given by density Pa{R) = {^\{R)'^ ct{,R)) ■ In accordance with fermionic 
commutation relations [15], the quasi-particle amplitudes are normalized as: 

j dR\uj{R)\^ + \vj{R)\'^ = 1 (4) 

and are related to the spin densities through : 

oo 

p^iR)=J2\^,iR)\'f{E,) 
i=i 

oo 

pm = J2\''j{R)\'fiEj), (5) 

where f{E) = l/[e^^^^^ + 1) is the Fermi-Dirac distribution function. In the unitary 
limit the Hartree terms gpa on the diagonal of Eq. ([2]) do not actually diverge but are 
unitarity limited |26]. How to incorporating the Hartree term in the unitarity limit 
is beyond the mean-field BdG formalism. For this reason we ignore this term in our 
calculations. The paring field or gap paramter is give by 

A(^) = g{i^^{R)^^{R)) = gJ2^,iR)v;{R)f{Ej) (6) 

j 

Since the Hartree terms are ignored in our analysis, Eqs. ([2]), (j4]) and constitute a 
closed set of nonlinear equations which we solve self-consitently. However as presented 
above, our formulation has one problem which arises from a nasty side effect of 
the contact interaction. The contact interaction assumes wrongly that all states are 
scattered in the same way regardless of their incoming energy and consequently sums 
in contributions from collisions at arbitrarily high energy which creates an ultra-violet 
divergence. Hence the gap equation [6] needs to be properly regularized as we now discuss. 



2.2. Regularizing the BdG equations 

Due to the assumption of contact interaction, the gap A is a function of the center-of- 
mass coordinate of the pair, R. To discuss the regularization, it is more convenient to 
re-introduce back the relative coordinate r, with which the gap is defined as 

A{R, f) = {^^{R + f/2)^^{R - f/2)) 

which diverges as when r — )• [24j. To regularize Eq. ([6]), one simply subtracts off 
the 1/r divergence to obtain the regularized equation [23]: 
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In practice, the convergene of the sum above is quite slow and we discuss here a 
numerically efficient way of evaluating A(i?) to sufficient accuracy without undue effort. 
First an energy cutoff Ec is used to break the sum of Eq. ([6]) into two pieces as ff. 

^— = 2^ u,{R)v^{R)f {Ej) + — ^ (8) 

The second term Ac(-R) is an approximation to the sum above the cutoff using the LDA 
result for the pairing field [2S] which can also be written as : 

A,{R) _ A{R) p d^k 

where is the momentum cutoff related to Ec- This leads to a computational efficient 
form of the gap equation: 

J2 u,{R)v]{R)f{E,). (10) 

UeE[,R) E,<Ec 



9 (27^)^ A, . [7^ 



Here we have employed the identity: 

dk e'^'' 1 



:iii 



/o (27r)3 1^ 27rr 

to subsume the LDA approximation of the gap (Ac) into an effective interaction defined 
by: 



1 



U,s{R) 



' d'k' 




2 T 



(12) 



Below this cutoff Ec, the quasiparticle states are calculated exactly by solving Eqs. ()2]), 
( fTOj) and ( fT2|) self-consistently along with the normalization conditions: 

N^= [ Pa{R)dR (13) 



which conserve total particle number = A^ + A| and overall polarization P = 
{N^ — Ni)/N. The iterative solution of these equations is achieved using a modified 
Broyden's approach [27] which is a nonlinear mixing scheme allowing the formation of 
polarized regions even if they were not present in the initial condition. In this scheme 
convergence was achieved when the root mean squared difference between A at different 

iterations was below some tolerance i.e., W ^ ^a^'^'"*"^ ^ ^^^^^ j position 
index and i represents the iteration number. We should point out that at unitarity this 
method is analogous to a descent technique where the step is optimized through the 
residuals stored from a few previous steps. As mentioned previouly we choose as our 
initial condition Alda, the LDA solution to the BdG equations. This is an important 
point becasue we found in previous work pi] that at larger particle numbers the solution 
of the BdG can be quite sensitive to the initial conditions. 
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2. 3. Special features 

We discretize using a linear triangular finite element mesh in the p-z plane which 
anticipates that our results will retain the cylindrical symmetry of the confining potential 
Vext- The accuracy of these calculations are controlled by the density of the trianglular 
mesh and the cut-off Ec used in the hybrid scheme. Both of these are changed in 
successive solutions until the free-energy or relevant observable converges to a sufficient 
accuracy. Experience has taught us that this simple renormalization scheme typically 
converges when the cutoff is of the order QEp (where Ep is the Fermi energy) which 
implies that the number of quasiparticle states to be directly calculated by Eq. (|2]) is 
about 6A^. Note that this puts a constraint on the density of the discretizing mesh. 
Thus, for moderate system sizes, we are still presented with a very large problem. For 
example, for N ~ 10^ particles, one essentially needs to calculate ~ 10^ quasiparticle 
states at each iteration. 

One important consequence of our finite element discretization is that it yields 
sparse matrices which are suitable to massively parallel matrix computations. This is 
of key importance given that the slow convergence of the sum in Eq. f lTOj) condemns 
us to calculate a very large number of quasi-particle states. This is true inspite of 
our efficient hybrid scheme, without which calculation would be prohibitive. It is 
immediately obvious that these difficulties will increase with the number of particles 
A^, and will make the problem impractical for even moderate particle numbers without 
very careful formulation. In our case these difficulties are inescapable since the issues 
to be addressed occur in the presence of finite size effects and confinement. Hence it 
was crucial to develop the ability to perform calculations with realistic particle numbers 
because it is not a priori obvious how physical properties will scale with system size. 
At each iteration, we need to find a large number of eigenvalues and eigenfunctions for 
large matrices. To this end, we use a novel shift-and-invert scheme which we developed 
independently but is very similar to the one described in Ref. |28j . 

Briefly, the scheme involves partitioning the sought spectrum amongst groups of 
processors working independently. The size of the group is determined by the minimum 
number of processors with enough total memory to store the inverted matrix which is 
required for building the local Krylov basis. The main challenge here is bookkeeping to 
prevent over-counting of states and balancing the load amongst the processor groups. It 
is conceivable that this method could have issues in cases where the equations support 
huge degenerate subspaces. In our particular formulation we exploited cylindrical 
symmetry and parity along the long trap axis to reduce the problem. Consequently 
we only had to contend with accidental degeneracies. A good analysis of these issues 
may be found in Ref. [28] but a more thorough description of the numerical details 
will be presented elsewhere. However we note here that this parallelization scheme is 
very efficient on distriubted computing systems and scales easily to thousands of CPU's 
which is as high as we have tested. Potentially it can be used to study much larger 
systems than we have reported here or in \21\ . 



Bogoliuhov-de Gennes study of trapped spin-imbalanced unitary Fermi gases 7 
2.4. Validity of the BdG 

We should devote a few lines here to comment on the validity of the BdG theory 
at unitarity becasue it is expected that quantum fluctuations and other effects due 
to the strong interactions could be significant in this regime. The main drawbacks 
of the BdG is that it fails to account for phase fluctations. At unitarity it has an 
additional disadvantage in that it also fails to account for interactions within the normal 
fluid which is unitarily limited [26]. However the BdG is widely expected to yield 
qualitatively reliable answers for two main reasons. First because of the finite size of 
these experiments, the trapped gas enjoys protection from fluctuations of arbitrarily 
low energy or of very long wavelengths. Secondly due to experimental evidence for 
superfluidity at unitarity, it is quite clear that interactions within the normal fluid are 
not so great that the order parameter cannot form or will be destroyed. Thus the 
failure to account for these effects is not expected to change the topology of the phase 
diagram but at most would slightly shift the phase boundaries. Since our purpose to 
examine the suitability of the LDA is qualitative, we are confident that the BdG can 
account for the essential physics. Nevertheless, due to the limitations within the BdG 
formalism, in particular the neglect of interaction in the normal fluid, our calculation 
fails to quantitatively locate position of the Clogston limit. In addition, it cannot be 
applied to study a system with extremely large polarization (i.e., P ^ 1) where the 
polaron physics will dominate [131 129] . 

3. Results for = 200 

In this section, we focus on a relatively small particle number of A^ = 200. As we 
will show, the system is rather sensitive to the trap geometry. In the following, we first 
briefly discuss the case of a spherical trap with a = 1 and then concentrate on elongated 
cigar-like trapping potentials with a > 1 and then discuss them in detail. 

3.1. Spherical trap 

Liu et al. studied a spin-imbalanced Fermi gas confined in a spherical harmonic trap in 
Ref. [16]. To benchmark our work, we first did a series of calculations for this geometry 
and found our results in perfect agreement with those reported in Ref. [16]. In this case, 
even though we anticipate only cylindrical geometry, the density profiles always obey the 
spherical symmetry. Note that the authors of Ref. [16] solved the one- dimensional (ID) 
radial equation, hence the spherical symmetry of the cloud is automatically imposed. 
We refer the readers to Ref. [16] for details; here we give just a brief description of the 
key features. The density profiles indicate a phase-separation scenario: a fully paired 
BCS superfluid core at the trap center surrounded by a fully polarized shell composed 
of excess majority spins. A thin layer of partially polarized gas forms the interface 
between the superfluid core and the normal shell. In this intermediate regime, the 
minority density and the order parameter sharply drop to zero. Here and in other cases. 
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Figure 1. Density profile of tlie atomic clound for N = 200 in an elongated 
trap with aspect ratio a = 5. The upper, middle and lower row are results from 
different polarizations P = 0.2, 0.4 and 0.6, respectively. In each row, we have 
shown (from left to right) the column densities of the majority component J dx pf, 
the minority component J dx pi, their difference J dx{p^ — pi), and the axial spin 
density nioiz) = J dxdy {pf - p^). 



we always found that the profile of the order parameter closely follows that of the density 
of the minority spin component. Furthermore, in this case, the LDA gives very good 
agreement with the full BdG calculation even for particle numbers as small as a few 
hundred. 

3.2. Cigar trap 

So far, all the experiments on spin-imbalanced Fermi gases have been performed in 
cigar-like traps with aspect ratio a > 1. For a given particle number, the ID regime 
is eventually encountered in this geometry by increasing a and creates the possibility 
to study the 3D-1D dimensional crossover. Figure [T] illustrates several examples of the 
density profiles for = 200 atoms confined in a moderately elongated trap with a = 5 
(this trap aspect ratio is close to what has been used in the MIT experiments). We find 
it convenient to express our results in terms of the Fermi energy Ep = (3iV)^/^a^/^, 
central number density (2£'^)^/^/(67r^), and the Thomas- Fermi radius along the ^-axis 
Zp = ^lEp for a single species ideal Fermi gas of N j 2 particles in a trap with identical 
parameters. The upper row of Fig. [1] shows the density profiles of a system with a 
relatively small polarization P = 0.2. Here the axial spin density nij:)[z) exhibits a 
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(a) P=0.2 




(b) P=0.6 



n- BdG 




Figure 2. Density and order parameter profiles along the axial and radial axes in a 
cigar-like trap with a = 5 for two differnet polarizations: (a) P — 0.2 and (b) P — 0.6. 



double- horn structure and vanishes near z = 0. This is a clear violation of the LDA 
which predicts that nic should be fiat topped [30]. Fig. [T]can be examined in tandem 
with Fig. [2] where for a closer inspection, we plot the densities and the order parameter 
along the axial and radial axis for two different polarizations. Fig. |2]^a) displays results 
for P = 0.2. The density profiles along the z-axis show clearly a phase separated three- 
region structure — moving from the center to the edge of the trap, we encounter a 
fully paired superfiuid core, a partially paired intermediate region and a fully polarized 
normal gas, just like in the previous case of spherical trap. In stark contrast, the density 
profiles for the two components along the p-axis are completely overlapped. In fact, this 
matching of the radial profiles occur for \z\ < 0.1. As a consequence, the axial spin 
density vanishes near 2; = as shown in the upper row of Fig. [H 

That the majority and minority densities overlap along the radial direction can 
be understood from an argument invoking the surface energy. When induced phase 
separation occurs, there is an accompanying surface energy associated with the interface 
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P-0.2 



P-0.4 



P = 0.6 




Figure 3. Same as in Fig. [U but for a = 50. 



between the two phases. The system will then try to minimize the interface in order 
to reduce the associated energy. For a cigar-like trap as we study here, the superfluid- 
normal gas interface area can be efficiently reduced if the two spin components match 
their densities radially. The authors of Ref. [311 [32] devised phenomenological theories 
to include the surface term variationally to explain the breakdown of the LDA observed 
in the Rice experiment O [7]. In our calculation, the surface energy is automatically 
included from the self-consistent BdG formulation |33l |. 

As polarization increases, eventually it becomes energetically unfavorable to have 
this radial overlap. This is illustrated in Fig. Mjo) for P = 0.6. Consequently, the axial 
spin density no longer vanishes near z = and the LDA becomes more accurate (see the 
middle and bottom rows of Fig. [1]). In addition, it is quite noticeable that, particularly 
for large P, the minority component density has a steeper down turn along the axial 
axis than along the radial axis. Moreover, in the partially polarized intermediate region, 
the order parameter has a small oscillation along the axial axis, but not along the radial 
axis. Similar order parameter oscillations were also found in the spherical trap case |16j . 
This is a consequence of the proximity effect which, in the context of superconductor, 
occurs when a superconductor is in contact with a normal metal, the Cooper pairs from 
the superconductor diffuse into the normal component. 

Next, we keep N fixed at 200 but increase the trap aspect ratio to a = 50 
which represents a much more elongated cigar trap and close to what is used in the 



Bogoliubov-de Gennes study of trapped spin-imbalanced unitary Fermi gases 



11 




Figure 4. a = 50 for two different polarizations: (a) P = 0.2 and (b) P = 0.7. Same 
units as in Fig. [2l We do not show the radial density and order parameter profiles, 
which look more or less like those in Fig. [2][a). (c) The surface plot of the order 
parameter for the case of P = 0.7 in the p-z plane, showing strong oscillations with 
the nodes aligned along the radial direction, (d) The density oscillations for P = 0.7 
leave a strong signal in the doubly integrated axial spin density uid- 

Rice experiment. A similar display of the column and axial spin density profiles for 
different polarization as in Fig. [T] is shown in Fig. [3l In this very elongated trap, the 
majority and minority components have their densities matched along the radial axis 
up to the highest polarization we have calculated which is P = 0.7, and the minority 
component has a boxy-looking density profile. This further confirms that the system 
is able to greatly reduce the effective surface area between the normal state and the 
superfluid state in anisotropic cigar-like traps. Another marked feature for such an 
elongated trap is the prominent oscillations of the order parameter along the z-axis. 
As demonstrated in Fig. Hj these oscillations are quite generic features in such a trap 
with finite P. As P increases, both the amplitude and the spatial extension of the 
oscillations increase. As shown in Fig. Hl^b), at large polarizations, the axial length 
of the partially polarized intermediate region becomes comparable to or even larger 
than that of the BCS core. Accompanied by the oscillation in the order parameter, 
the density profiles (in particular, the minority density) also exhibit strong oscillations. 
Such oscillations are reminiscent of the FFLO pairing state predicted by Fulde, Ferrel, 
Larkin and Ovchinnikov [Sj H] in which the Cooper pairs possess finite momentum and 
the order parameter in the bulk develops sinusoidal oscillations that break the spatial 
translation symmetry. 

Our calculations also show that these axial oscillations are aligned along the radial 
shown in Fig. Hl^c). We even intentionally started from an initial ansatz of 
A where the axial oscillations are present but with the nodes mis-aligned in the radial 
direction, the BdG iterations eventually converge to a state where the nodes are perfectly 
aligned radially. This radial alignment has important impact in detecting the oscillations 
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in column density profiles where the densities are integrated along one radial axis: Due 
to the radial alignment, the oscillations are not washed out and can be easily observed, 
for example, in the doubly integrated axial spin density nm, as illustrated in Fig. Hl^d). 

It is interesting to compare our result with the recent work by Bulgac and Forbes 
[23] who, using a density functional theory, argued that the FFLO pairing phase occupies 
a larger phase space region than people previously thought for a 3D homogeneous 
system. The FFLO state found in Ref. |23] is also associated with large-amplitude 
density oscillations, particularly in the minority component. Another perspective on 
the order parameter oscillations and its potential connection with the FFLO phase is 
dimensionality. It is well known that the partially polarized phase with FFLO-like 
oscillations is prominently featured in the phase space of ID systems [25| IMf 135] . That 
the reduced dimensionality favors such a state can be understood from the argument of 
Fermi surface nesting or alternatively from the cost of creating domain walls. The cigar- 
like traps used in our calculation mimic a quasi- ID system and may be the reason that 
we see pronounced oscillations in our calculation. If this latter explanation is correct, 
i.e., the partially polarized region featuring FFLO-like oscillations is due to the effective 
reduction of the spatial dimension, we then expect to see these oscillations diminish as 
is increased while the trap aspect ratio is fixed, which makes the system more 3D-like. 
To confirm this, we now turn to the next section where we keep a = 50 but vary the 
total particle number N. 

4. Results for large particle numbers at a = 50 

We now consider a trapped system in a very elongated cigar trap with a = 50. Fig. [5] 
shows the density and order parameter along the axial and radial axis at a fixed 
polarization P = 0.3 but different values of total particle number A^. As one can clearly 
see, the oscillations in both the order parameter and the density profiles diminish as 
A^ is increased and the LDA approximation becomes more and more accurate, which 
indicates that the FFLO-like region observed above for small A^ does not represent a 
bulk 3D phase. Rather, it is a finite-size effect due to the effective reduction of the 
spatial dimension. 

Nevertheless, we have discovered that as A^ increases, the system exhibits a tendency 
towards metastability [21]. Numerically, by starting from different initial anstze for the 
order parameter A(r^, the BdG solution may converge to different final states. Our 
calculations show that among these different states, the one that closely resembles the 
LDA solution always has the lowest energy as long as A^ is sufficiently large (A^ > 10"^), 
but that there may exist several metastable states with energies just slightly larger 
that violate the LDA. The observed LDA-violating states at Rice are most likely these 
metastable states. Experimentally, whether the ground state or a metastable state will 
be realized may depend upon how the evaporative cooling procedure is implemented 
[36] . This has been confirmed very recently in a new experiment by Hulet group [37] . 
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(a)A^= 1100 




ih)N= 10000 
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(c) iV= 50000 
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Figure 5. Density and order parameter profiles along the axial and the radial axis in 
a cigar-like trap with a = 50 for P = 0.3 but different values of total particle number 
N. Same units as in Fig. [21 
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5. Conclusion 

In conclusion, we have carried out a systematic study of a trapped spin-imbalanced 
Fermi gas in the unitary limit up to a total number ~ 10^ atoms. We study a 
class of solutions which has recently been identified as having the lowest energy [21], [22] 
through a self-consistently solution of the Bogoliubov-de Gennes equations using state- 
of-the-art numerical techniques. For a given set of trapping parameters, the LDA will 
eventually become accurate for sufficiently large A^. However, the validity of the LDA 
is also sensitive to the trap geometry: Traps with small anisotropy favor the LDA. 
Our calculations show that for a relatively small number of atoms in a very elongated 
trapping potential, the system contains three phases: an unpolarized BCS phase, a 
partially polarized FFLO-like phase and a normal phase. That the FFLO region exists 
may be understood from the view point of reduced effective spatial dimension. As A^ 
is increased while all other parameters remain fixed, the FFLO-like region eventually 
disappears. A detailed analysis of dimensional crossover will be reported in the future. 
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